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This paper is devoted to the study of the dynamics of two weakly-coupled Bose-Einstein con- 
densates confined in a double- well trap and perturbed by random external forces. Energy diffusion 
due to random forcing allows the system to visit symmetry-breaking states when the number of 
atoms exceeds a threshold value. The energy distribution evolves to a stationary distribution which 
depends on the initial state of the condensate only through the total number of atoms. This loss of 
memory of the initial conditions allows a simple and complete description of the stationary dynamics 
of the condensate which randomly visits symmetric and symmetry-breaking states. 
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I. INTRODUCTION 

The experimental achievements of Bose-Einstein condensates (BECs) with various dilute trapped alkali-metal atom 
gases have opened new possibilities to investigate macroscopic quantum effects lj. Indeed it is well known that 
phase coherence plays a crucial role in macroscopic quantum systems, and different systems have been addressed: 
Josephson effect in superconductors Q, weakly linked superfluid He-B reservoirs two coupled BECs in a double 
well trap Q, and coupled condensates in different hyperflne levels In this paper we shall address the quantum 
coherent tunneling between two BECs trapped in a double-well potential, and we shall focus our attention to the 
influence of random external fluctuations. In superconductors the main origin of external noise is thermal or quantum 
fluctuations. Mathematically the problem can be reduced to the pendulum equation (sine-Gordon equation) with 
external fluctuations <fitt + sin(0) + 70t = f(t) [|| and the energy diffusion of the system can be studied. In particular 
Stratonovich has studied this problem in the strong damping case jcfit + sin(0) = f(t) [3. 

In BECs new types of problems appear due to the existence of the trap potential. Coherent phase phenomena in 
coupled BECs are described by a system of ordinary differential equations for the imbalance in atomic population z 
and relative phase <fi. This system can be considered as the equation for the nonrigid pendulum. The pendulum case 
recovered if the imbalance is small z 2 « 1. Nonrigidity leads to new phenomena, and we can distinguish two regimes 
in the atomic imbalance evolution. The first regime is the macroscopic quantum tunneling (MQT), when the atomic 
population is periodically exchanged between wells with (z) = 0, which is analogous to the Josephson effect. The 
second regime is the macroscopic quantum localization (MQL) phenomena, when the mean field nonlinearity excesses 
some critical value, the mean value of the imbalance becomes nonzero (z) ^ 0, which corresponds to the localization 
of the atomic population in one of the wells in the form of small-amplitude oscillations around the bottom of the well. 

It is therefore interesting to investigate the influence of fluctuations on the evolution of the system in these states. 
In particular we would like to pay a particular attention to the possible switching between these states. This problem 
has a general character and it is important for the theory of nonlinear directional couplers in nonlinear optics ja, and 
nonlinear dimers |icj|. Periodic variations of the parameters of the trap have already been considered, and resonant 
phenomena have been exhibited for the weakly-coupled BEC an d for the strongly overlapped BEC The 
influence of a periodic time- varying atomic scattering length was addressed in [l3j |. Dynamical tunneling between 
regions of regular motion was shown to be involved by strong periodic modulation of the tunnel coupling between 
the two modes jl4| . Finally macroscopic quantum chaos driven by a time-periodic trap asymmetry was predicted in 
|15| . In this work we shall study the influence of random fluctuations on the BEC dynamics trapped in a double- well 
trap. We shall address the case of fluctuating zero point energies so that the energy difference AE(t) is a zero-mean 
random process. Such variations are induced by small oscillations in the barrier laser position |l6l Il7| . Indeed, the 
barrier is generated by a Gaussian laser sheet focused near the center of the harmonic trap and small oscillations of 
the barrier-laser position introduce the fluctuations of the zero-point energies. 
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II. BOSE-EINSTEIN CONDENSATES IN A DOUBLE- WELL TRAP 

The problem of BECs in a double-well time-dependent trap is described by the Gross-Pitaevskii (GP) equation 

ih% = -^-A* + V tr (r,t)* + g\*\ 2 V, (1) 
2m 

where Vt r is the double-well potential, g = 4irh 2 a s /m, and a s is the atomic scattering length. For weakly overlapped 
condensates (high barrier or well separated wells) we can use the two-modes decomposition 

* = ^i(t)*i(r)+^(t)*a(r), (2) 

where ^1,2 are complex time-dependent amplitudes of condensates in wells and are approximate ground state 
solutions of GP equation in first and second wells respectively. Substituting this solution into ||TJ, multiplying the 
equation by $i 9 and integrating over the spatial variable we obtain the system of equations for the two modes 

M*),Mt) mm 



iK^ = [E l {t) + ai\^ l \ 2 ]i> l -K'ip 2 , (3) 
at 



ih^ = [E 2 {t)+aM 2 ]i>2-K'ip l , (4) 
at 

where E\^ 2 are the zero point energies, K is the coupling constant between the two modes, and a\ t % are the mean 
field nonlinearities for both modes. Here we take into account a small oscillation in the laser-barrier position so that 
AE = Ei — E 2 is time-varying [2(J. We write tpj — y/Nj exp(idj). N T — Ni + N 2 is the (constant) total number of 
atoms. The fractional population imbalance 

. « (5) 

and the relative phase 

<f>(t) = 9 2 (t) - e^t) (6) 

satisfy 



zt = - \J\ - z 2 sm{4>) - i)4>u (7) 

fa = -AE(t) + Az + JL— Cosi®, (8) 
V 1 — z z 

where we have rescaled to a dimensionless time t2K/h — > t and we have introduced new variables 

A F ffl _ E^-E^t) OL\ OL 2 

AE{t) 2K + ~^k~ NT ' (9) 

aN T a\ +a 2 , , 

A = ^T' a = —2— (10) 

Note that we have included the damping term —r\<\>t in Eq. (7J that takes into account a non-coherent dissipative 
current of normal-state atoms, proportional to the chemical potential difference 

III. HAMILTONIAN STRUCTURE OF THE UNPERTURBED SYSTEM 

The Hamiltonian of the unperturbed system (AE = 0, r\ = 0) is 



E = — - xjl-z 2 cos(4>). (11) 

E can take any value between E m i n = — 1 and E max = 5(A+ 4). It is an integral of motion. The orbits of the motion 
are closed, corresponding to periodic oscillations. In order to explicit the periodic structure of the variables z and <fi, 
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we introduce the action-angle variables. The orbits are completely determined by the energy E imposed by the initial 
conditions. The period is denoted by T(E), the motion is described by 

z(t) = Z(E,0(t)), (12) 
cos(cl>(t))=C(E,6(t)), (13) 
sra(0(t)) =S(E,6(t)), (14) 

where Z, C, and S are smooth functions, periodic with respect to 9 with period 2n, the angle satisfies 

dO 2tt 



dt T(E) ' 

and the action is defined by 



I{E) = h fj( x ) dx - 



(15) 



(16) 



The analysis is standard but quite cumbersome, so we only list the final results that can be written in terms of 
elliptic functions. We need to introduce a series of parameters. 

I p_1 A 2 I 

k(E) = ——, aE) = — r~ k{E) (17) 



M(E) = ~(l + ^y z 2 (E) = jyJ K +y/C. (18) 



Three cases can be distinguished. 

1. If E £ [-1, 1), then M < 1, the period is 



T(E) = ^K(M), (19) 



where K is the complete elliptic function (22 



c( "> = iff=f' (21) 

M0 -. J^ -d.fHWg,^. (22) 



2-/MV1 - Z 2 

cn, sn, and dn are tabulated Jacobian functions j^. Note that Z is an even function with respect to 6. These solutions 
preserve the z-symmetry in the sense that they correspond to eigenf unctions of the Gross-Pitaevskii equation that 
are either odd or even functions. The stationary ground state has energy —1 and it is given by z = and <j) = 0. For 
E £ (—1, 1), we have (z) — and the atomic population is periodically exchanged between the two modes. This is 
the MQT regime. 

2. If E £ (l,E max ), then M > 1, the period is 

T(E) = £-K(±). (23) 
Az2 M 

The three periodic functions Z, C, and S are given by 

Z(E,0)^±z 2d n^) 6 -^), (24) 
K? 2 — 9E 

C{E,8) = ±±-^Z, (25) 
V ' 2^/Y^Z 2 

"z 2 ( „, 1 N 6 1 



S(E,9) = sncn K(— )-,— . (26) 
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FIG. 1: Damping term E Af(E) for A = 2 (picture a) and A = 10 (picture b). 



Note that Z is either always positive, or always negative-valued. These solutions correspond to eigenfunctions that 
break the z-symmetry. In particular there exist two stationary solutions with the maximal energy E max that are 
z = ±-v/l — 1/A 2 , <f> = 7T. Their existences result from the nonlinear interatomic interaction and it is possible only 
when the number of atoms is large enough so that nonlinear coefficient A > 1. If E e (l,E max ) the population 
imbalance periodically oscillates around a non-zero average value (z) ^ 0. This is the MQL regime. 

3. If E = 1, then M = 1, the motion is given, in this special case, by the non-oscillatory hyperbolic secant 

z{t) = z 2 sech (yi - l(t + t )) , (27) 
with z 2 = 2\/A - 1/A and t = (l/VA— l)argch(z 2 /z ). 



IV. EFFECTIVE DYNAMICS 



A. Damping 

In presence of a small damping an adiabatic approach is possible which yields that the energy decays as 

dE 

— = -vN{E), (28) 

where 

N{E) = 1- ( AZ + ■ z 

AT(E) vanishes at E = —1, E = 1, and E = E max . It takes large values just above and below E = A/2. In Figure^] 
we plot the function E Af(E) for two different values of A. Expansions of AT are also presented in the Appendix. 

Damping is not able by itself to break the symmetry as it involves an energy decay, while symmetry breaking 
occurs when the energy goes from a value below the separatrix 1 to a value above the separatrix. If the system 
starts from one of the two breaking states with the energy E max then it stays exactly at this state as N(E max ) = 0. 
However, a linear stability analysis shows that these states are not stable. More exactly, if the initial state has 
energy E max — e with e«l, then it stays close to this state during a time of the order of ~ [A(A 2 — l)e] _1 because 
N(E max — e) ~ A(A 2 — l)e; after this time it quits the breaking state and converges to the ground state E m i n . This 
behavior is described in Figure Note the remarkable agreement between the results from numerical simulations of 
the system of equations Ij7l8(l with AE = and the effective solution given by Eq. (J2EJ)- 



(E, 0)dO. 



(29) 
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FIG. 2: Picture a: Evolution of the population imbalance z in absence (dashed line) and in presence (solid line) of a small 
damping r\ = 10~ 2 . Here A = 2 and the initial state is zo = v/1 — 1/A 2 — 0.02 ~ 0.846 and 0o = t- Picture b: Evolution of 
the energy E(t) = \z 2 (t)/2 — yl — z 2 (£) cos(0(t)) according to the numerical resolution of the system I|7I8[I (solid line) and 
according to the effective equation 1281 1 (dashed line). 



Random fluctuations 



Let us consider a random perturbation of the system that can be written in the form 



Zt 



d z H, 



(30) 
(31) 



where H = Ho(z,<f>) + m(t)V(z,cos(<f>),sm(<f))), Ho is the unperturbed Hamiltonian Qllfl. and m(t)V is the time- 
dependent perturbation, m is assumed to be a stationary, zero-mean, random process. We think in particular at a 
time-dependent AE(t) in Eq. Q so that m(i) = AE(t) and V = —z. 

In presence of perturbations, the motion of (z, <f>) is not purely oscillatory, because the energy and the action are 
slowly varying in time. We adopt the action-angle formalism, because it allows us to separate the fast scale of the 
locally periodic motion and the slow scale of the evolution of the action. The motion is governed by the system 



It = -m(t)hg(I,6), 

9 t = Lj(I) + m(t)hi(I,9), 



(32) 
(33) 



where 



2tt 



I i ► £(I) is the inverse function of E ^I(E), and h(I,9) = V(Z(I,9),C(I,0),S{I,0)). Using 



To£(I) ' 

standard diffusion-approximation theory |21|. we get that the action / behaves like a diffusion Markov process with 
the infinitesimal generator 



c 1 d 

Cl -2dl 



The diffusion coefficient is 



MX) = ~ / Mi", 0)hg{1, 9 + uj(I)t) (m(O)m(t)) dtd9, 

7T Jo JO 



(34) 



(35) 



where the brackets stand for a statistical averaging. This means in particular that the probability density function of 
I(t) satisfies the Fokker-Planck equation dtp = C}p, p(t — 0, 1) — 5(1 — Iq), where Iq is the initial action at time 
and Cj is the adjoint operator of Ci, which is equal to Ci in our configuration as Ci is self-adjoint. 



C. Energy diffusion 



The results of the two previous subsections can be combined to address the case of system H7I8JI with a white 
noise model for AE and a damping rj > 0. We get that the energy E of the system is a diffusion process with the 
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FIG. 3: Period T(E), drift and diffusion coefficients for A = fO. 



infinitesimal generator 



where 



2TT 2 a d 
T{E)dE 



B(E) 8 
T{E)dE 



B(E) 



1 

2^ 



a 



(AE(O)AE(t)) dt. 



'-de, 



(36) 

(37) 
(38) 



The infinitesimal generator is not self-adjoint. The drift of the energy is d(E) = a2n 2 d E [B(E) /T{E)]/T{E) - r]Af(E) 
and the diffusion coefficient is a(E) = a2ir 2 B (E)/T 2 (E). They are plotted in FigureOHm the case 77 = 0. This result 
allows us to compute all relevant quantities, in particular the stationary energy distribution. Indeed the diffusion 
process E(t) is ergodic, and the statistical distribution of the energy becomes independent of the initial state for large 
time t. ft converges to a stationary distribution that can be computed explicitly as the solution of the elliptic equation 
£%P = 0. 



V. ENERGY DISTRIBUTION 

The stationary energy distribution for the system is determined by the resolution of the equation C* E p = that is 

rE T 2 N \ 



p(Et) = CT(E)exp(-l 1^ 



2ir 2 B 



(x)dx , (39) 



where the normalization constant C is chosen so that f E ™ ax p(E)dE = 1. This distribution depends only on A and 
the ratio rj/a. A first important quantity is the average energy 

(E) = J Ep(E)dE. (40) 

It is plotted in Figure^, as a function of A for different damping rates, which shows in particular that the average 
energy grows linearly with A when r\ = 0, according to E ~ A/6, but possesses a maximum if r\ > 0. Another 
important quantity is the proportion of time spent by the system in breaking states 

Rb = J p(E)dE, (41) 

which is zero for A < 1 and becomes positive for A > 1. Rf, is plotted in Figure^ which shows that, for a given 
value of 77/a, there exists a critical value for the nonlinear coefficient A such that the proportion of time spent by the 
system in breaking states is maximal. We can examine theoretically several particular cases. 
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FIG. 4: Picture a: Average energy as a function of A for different damping rates. Picture b: Proportion of time spent by the 
system in breaking states Rt as a function of A for different damping rates: rj/a = (thick solid line), -q/a — 1CP 5 (thin solid 
line), rj/oL = 10~ 4 (dashed line), r//a — 10 -3 (dash-dotted line), 77/0 = 10~ 2 (dotted line). 



0.04 




FIG. 5: Proportion of time spent by the system in breaking states Rb as a function of A. The numerical integration of Eq. I|41|l 
is plotted in solid line, the approximate formula Rb — 0.21A 3//2 is plotted in dashed lines. 



A. Absence of damping 



In absence of damping r/ — 0, the energy distribution does not depend on the diffusion rate and it is given by 

p[E ) = -^1^1 . (42) 

Let us first address the case where the number of atoms is just above the critical number that allows the existence of 
symmetry-breaking states. If A = 1 + A, Q < A <C 1, then we get by using the expansions described in Appendix IA 41 
that 



4( 



id)7^(5-^r)> * £ €[-1,1), 



PiE) ^ * r*fc<gra. if*€(i,^, (43) 



c[A+ A /2(B max -_E)] 1 /2 \X+y/2(E maa: -E) 

with c ~ 4.44 and the proportion of time spent by the system in breaking states is Rb — 0.21A 3 / 2 . A comparison 
with the numerical integration of 141|) shows that the approximate formula Rt, ~ 0.21A 3 / 2 is valid for < A < 0.2 (see 
Figure [5J. Thus, for a small number of atoms, it is not likely to observe the MQL regime. 

Let us now consider the case of a large number of atoms. If A ^> 1, then we get by using the expansions described 
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in Appendix I A 31 that 



^K^), if E& [-1,1 



m*l™i i 1 . E ( 44 ) 



V2A(S-1) : 



and the proportion of time spent by the system in breaking states is Rb — 1 — 2.55A -1 / 2 . We can see in Figure^ that 
this approximate formula is valid for A > 40. Thus Rb is close to one which means that it is very likely to observe the 
MQT regime. We can also compute the transition time from one breaking state to the other one. This computation 
is based on a well-known result of stochastic analysis which claims that fiEo.EtiE), the mean time to exit the interval 
[Eo, Ei] starting from the energy E, satisfies Ce^ = — 1 with the boundary conditions h(Eq) — (a(Ei) — 0. The result 
is that the transition time is of order A 2 /a. The computation of the transition time is similar to the one performed 
to obtain the solution to the Kramers' exit problem which is concerned with noise activated escape from a potential 
well |23j. As a result we can predict that, for large A, the system randomly goes from one breaking state to the other 
one according to a Markovian dynamics with an average period of A 2 /a. 

B. Small damping 

In presence of damping Eq. I|39|l shows that the stationary energy distribution is shifted toward the low energy 
region. To illustrate this phenomenon, we can consider and analyze the case A > 1 and rj/a -C 1/A. Using the 
expansions described in Appendix I A 31 the energy distribution can be written as 

2£^K(^), if£e[-i,i), 
,<£ll %M- ! ?4 (45) 



with 




if r]/a < 1/A 3 
f^^T if 1/A 3 « n/a « 1/A ™ 



As a result, the proportion of time spent by the system in breaking states is Rb — 1 — 2.55C, ; A 1 ' 2 which decays with 
rj. It is less and less likely to observe the MQT regime as the damping becomes larger. 

C. Strong damping 

Assume that rj/a ^S> 1/(1 + A). Using the expansions derived in Subsection IA li the diffusion operator can be written 



as 



^ = ^(1 + £)^(1 + A)(1 + £)A. (47 ) 



The stationary energy distribution is 

p(E) = i-exp (-^r^) M-hE max ]( E )> (48) 

where E c = a/[(l + A)t/]. The average energy is —1 + E c . The proportion of time spent in breaking states is 
exponentially small Rb — exp(—2/E c ). 

D. Numerical experiments 

In this section we compare our theoretical predictions with numerical simulations of Eqs. I)7I8|I . The random 
fluctuations of AE(t) are modeled by a stepwise constant process 



AE(t)=a'£x j l UtctU+1)ta) (t), 
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FIG. 6: Picture a: A typical z-trajectory for the set of parameters (a, t c ) — (4,0.05). Picture b: Average energy versus time 
for three different configurations: (<r,t c ) = (4,0.05) (solid line), (o,t c ) — (3,0.1) (dotted line), and (a, t c ) = (2,0.05) (dashed 
line). 

where the Xj are independent and identically distributed random variables with uniform distribution over (— </3, \/3) 
and t c is the coherence time. The coefficient a is then given by 

a = . 

2 

Damping is absent r\ = and the parameter A = 15. Three series of simulations are performed with the parameters 
(cr,t c ) = (4,0.05), (cr, t c ) = (3,0.1), and (a, t c ) = (2,0.05). Note that the first two configurations give almost the 
same value for the effective parameter a ~ 0.4. We have carried 1000 simulations for each configuration. The initial 
conditions are z$ = 0, <f>o = 0, but we have checked that the initial conditions play no role in the long-time dynamics. 

In Figure^ we plot a typical z-trajectory for the configuration (a, t c ) — (4, 0.05). As predicted by the theory, noise 
induces an energy diffusion which makes the system visit symmetry-breaking states. More quantitatively, the system 
spends around half-time is symmetry-breaking states as predicted by Figure QJa, and the life-time of a symmetry- 
breaking state is of the order of A 2 j a ~ a few hundreds. 

In figure 03 we plot the energy of the system for the three configurations. The energy is averaged over the 1000 
simulations. We can observe that the behavior is the same for the first two configurations, which is in agreement with 
the fact that both configurations possess the same effective parameter a ~ 0.4. We can also compare the asymptotic 
value of the average energy with the theoretical one (E) ~ A/6 = 2.5. Note that the time necessary to reach the 
asymptotic value is of the order of A 2 /a ~ 500. The configuration (<r,t c ) = (2,0.05) corresponds to a smaller value 
of a = 0.1. As a result the time necessary to reach the asymptotic energy value is longer, but the asymptotic value 
itself is the same as for the other configurations as it depends only on A. 

VI. CONCLUSION 

In this paper we have considered the dynamics of two weakly-coupled BECs. We have investigated the influence 
of small random oscillations of the barrier-laser position which induce fluctuations of the zero-point energies of the 
double-well trap. We have shown that random noise induces an energy diffusion which allows the system to visit 
symmetry-breaking states when the number of atoms (proportional to A) exceeds a threshold value. We have shown 
that the energy distribution evolves to a stationary distribution which depends on the initial state of the BEC only 
through the total number of atoms, and not through the initial imbalance or phase difference. We have computed 
the time necessary for the loss of memory of the initial conditions and the establishment of the stationary dynamics. 
Then we have described the stationary dynamics of the BEC which visits symmetric and symmetry-breaking states 
according to a random Markovian dynamics. We have investigated the role of damping in this process. In particular 
we have shown that, for a given damping, there exists a critical value for the number of atoms where the proportion 
of time spent by the BEC in symmetry-breaking states (i.e. in a macroscopic quantum localization state) is maximal. 
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The minimal energy is E min 

M - 2Tiw e ' Z2 - 



APPENDIX A: EXPANSIONS 

1. E Close tO E min 

■1. We can expand all quantities when E = E min + e, e > 0. We get in particular 



so that 



T(E) = 
Z{E,6) 



2tt 



v/T + T 
V2 



V / ecos(0), 



\/TTA 
C(£,0) = l-esin(<9) 2 



Af(£) = (1+A)e, 



(Al) 

(A2) 
(A3) 

(A4) 
(A5) 



2. i? close to _B„ 



Here we assume that A > 1 so that E max > 1. We can expand all quantities when E = E max — e, e > 0. We get in 

2(A 2 -1 

27T 



particular M ~ z 2 ~ ^1 + , 

A ^/A(A 2 -l) ' 



T{E) 



VA^T' 
VP - ! 



7 2e 



so that 



ZIE, 6) = 

A ^A(A 2 - 1) 

C(£,0) = -1 + Aesin 2 (6»), 



N(E) = A(A 2 - l)e, 

s < £ > - Atitry 



cos(#), 



(A6) 

(A7) 
(A8) 

(A9) 
(A10) 



3. Large A 

Let us assume that A > 1. If E e [-1, 1), then M ~ (1 + £)/2 < 1, 
If £ e (l,E max \, then by setting £ = Ae, e e (0, 1/2], we have 

V^TT 



T(25) 



A^' 



•Z(£,6>) = V2^ 



^V^ Sm ( 2 } 



(All) 
(A12) 

(A13) 
(A14) 
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so that 



N{E) = 2A 2 e, (A15) 

*w = (A16) 



4. A just above 1 



< A « 1. If E e [-1, 1), then M ~ i(l - and 



This case is interesting in that it is the regime where symmetry-breaking appear. Let us assume A = 1 + A, 

_ 

2\/2 

Note that E max = 1 + A 2 /2 + 0(A 3 ). If £ = 1 + ^e, e e [0, 1], then M ~ ±(1 + ^/==) and 

T(£) = - 2 ^ gf^V e' = l-. 



T(£\) = ^Jf Q(l-V^)), £' = (l-£)/2. 
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